Phylogenomics of darkling beetles (Coleoptera: Tenebrionidae) from the Atacama Desert

Background Tenebrionidae (Insecta: Coleoptera) are a conspicuous component of desert fauna worldwide. In these ecosystems, they are significantly responsible for nutrient cycling and show remarkable morphological and physiological adaptations. Nevertheless, Tenebrionidae colonizing individual deserts have repeatedly emerged from different lineages. The goal of our study was to gain insights into the phylogenetic relationships of the tenebrionid genera from the Atacama Desert and how these taxa are related to the globally distributed Tenebrionidae. Methods We used newly generated transcriptome data (47 tribes, 7 of 11 subfamilies) that allowed for a comprehensive phylogenomic analysis of the tenebrionid fauna of this hyperarid desert and fills a gap in our knowledge of the highly diversified Tenebrionidae. We examined two independent data sets known to be suitable for phylogenomic reconstructions. One is based on 35 neuropeptide precursors, the other on 1,742 orthologous genes shared among Coleoptera. Results The majority of Atacama genera are placed into three groups, two of which belong to typical South American lineages within the Pimeliinae. While the data support the monophyly of the Physogasterini, Nycteliini and Scotobiini, this does not hold for the Atacama genera of Edrotini, Epitragini, Evaniosomini, Praociini, Stenosini, Thinobatini, and Trilobocarini. A suggested very close relationship of Psammetichus with the Mediterranean Leptoderis also could not be confirmed. We also provide hints regarding the phylogenetic relationships of the Caenocrypticini, which occur both in South America and southern Africa. Apart from the focus on the Tenebrionidae from the Atacama Desert, we found a striking synapomorphy grouping Alleculinae, Blaptinae, Diaperinae, Stenochinae, and several taxa of Tenebrioninae, but not Tenebrio and Tribolium. This character, an insertion in the myosuppressin gene, defines a higher-level monophyletic group within the Tenebrionidae. Conclusion Transcriptome data allow a comprehensive phylogenomic analysis of the tenebrionid fauna of the Atacama Desert, which represents one of the seven major endemic tribal areas in the world for Tenebrionidae. Most Atacama genera could be placed in three lineages typical of South America; monophyly is not supported for several tribes based on molecular data, suggesting that a detailed systematic revision of several groups is necessary.


INTRODUCTION
Tenebrionidae Latreille, 1802 (Insecta: Coleoptera) have a worldwide distribution and are one of the larger families with more than 30,000 described species (Bouchard et al., 2021). In the majority of species, both larvae and adults are detritivores and often play a significant role in terrestrial food webs (Matthews et al., 2010). Based on their ecological preferences the Tenebrionidae can be broadly divided into two groups: species associated with trees and species with a shift in larval habitat from decaying trees to soil (Matthews et al., 2010). The latter are widely recognized as the insect group best suited for colonizing arid environments and are found worldwide in desert ecosystems. They have developed numerous morphological, physiological and behavioural adaptations to cope with extremely arid conditions and are therefore largely responsible for most of the nutrient cycling in deserts (Cloudsley-Thompson & Chadwick, 1964;Cloudsley-Thompson, 2001;Crawford, 1982;Matthews, 2000;Matthews et al., 2010;Cheli, Bosco & Flores, 2022;Raś, Kamiński & Iwan, 2022). Different from most other insect groups, their biodiversity sometimes increases with aridity (Kergoat et al., 2014a;Koch, 1962;Pfeiffer & Bayannasan, 2012). The genetic basis for these desert adaptations is not yet clear, but it is known that different lineages of the Tenebrionidae have repeatedly migrated into developing deserts in a convergent scenario (Matthews et al., 2010). Currently, 11 subfamilies, 106 tribes and 2,307 genera of Tenebrionidae are recognized (Bouchard et al., 2021), mainly based on the morphological characters (Doyen, 1972;Doyen, 1993;Doyen & Tschinkel, 1982;Kamiński et al., 2020;Matthews et al., 2010;Watt, 1974).
Recent analyses in insect phylogeny resolved the higher-level relationships in many cases using extensive molecular datasets (e.g., Chesters, 2020; Misof et al., 2014;Wipfler et al., 2019). The intra-ordinal relationships in Coleoptera (Bocak et al., 2014;Cai et al., 2022;Gunter et al., 2014;Hunt et al., 2007;McKenna et al., 2019;Zhang et al., 2018) and the intra-familial relationships of the larger beetle families (e.g., Tarasov & Dimitrov, 2016;Nie et al., 2020;Shin et al., 2018;Souza et al., 2020) was also the focus of several such studies. Regarding the Tenebrionidae, unresolved relationships were repeatedly addressed by molecular analyses in recent years, which, among others, consistently confirmed the monophyly of the family (Gunter et al., 2014;Kergoat et al., 2014b;Kamiński et al., 2020). However, these phylogenetic reconstructions are still under discussion because the internal relationships are still not fully solved. In particular, the subfamilies Tenebrioninae Latreille, 1802 and Diaperinae Latreille, 1802 appear to be artificial groups that require thorough revaluation. (e.g., Aalbu et al., 2002;Kergoat et al., 2014b;Kamiński et al., 2020;Johnston et al., 2020). A recent study convincingly suggested the subfamily Blaptinae Leach, 1815 as a monophyletic group based on molecular and morphological analyses (Kamiński et al., 2020); this lineage contains taxa that have traditionally been placed within the presumably polyphyletic subfamily Tenebrioninae. One of the limitations of all of these phylogenetic reconstructions is the lack of comprehensive sampling of lineages from southern Africa and southern South America. Both Southern Africa and South America each have a highly conspicuous tenebrionid fauna including several endemic tribes (e.g., Carrara & Flores , 2015;Koch, 1962;Kuschel, 1969;Matthews et al., 2010;Kamiński et al., 2021) and contain two of the oldest and driest deserts in the world, the Namib and Atacama Deserts (Clarke, 2006;Goudie & Eckardt, 1999) where tenebrionids represent one of the most abundant insect groups.
Aridity in the Atacama Desert can be traced to the Triassic, but the current conditions are closely related to the Andes uplift in the Miocene (Clarke, 2006), because this mountain range acts as an effective rain shadow (Houston & Hartley, 2003). The regions west of the Andes experienced a long-term decrease in precipitation in this context; the corresponding aridification presumably started in the early Miocene in what is now the core area of the Atacama Desert (Dunai, González-López & Juez-Larre, 2005;Ritter et al., 2018) and intensified throughout the Miocene until the present (Jordan et al., 2014;Ritter et al., 2018). Today, the core of the Atacama Desert (Central Depression between 19 • S−23 • S) is characterized by hyperarid conditions with less than two mm/yr of precipitations (Houston, 2006), making it one of the driest regions on Earth (Clarke, 2006). These climatic conditions are apparently a barrier for the evolution of organisms, and even well-adapted xerophilous insects as darkling beetles avoid the core of the Atacama Desert. Indeed, most tenebrionids prefer peripherally located and slightly wetter habitats in the Coastal and Andean Cordilleras (Fig. 1). However, the long-lasting interactions between tectonic activity and past climate changes in the Atacama Desert created conditions for the diversification of a very peculiar fauna of tenebrionids, some with very ancient relationships (see Endrödy-Younga, 1996;Ferrer, 2015); and under the influence of the fauna of neighboring regions of the Peruvian Desert and the Intermediate Desert of Coquimbo (Peña, 1966a).
The main goals of the current study are obtaining insights (1) into the phylogenetic relationships of the Atacama genera and (2) of the relationships of these taxa to Tenebrionidae from other regions. For this purpose, we collected material for molecular analyses of almost all described tenebrionid genera (30 genera including an undescribed genus of Alleculinae Laporte, 1840) that inhabit the Chilean Atacama Desert including the adjacent Andean Cordillera. Since it is unlikely that analyses of individual genes can resolve all issues concerning the higher phylogeny of the Tenebrionidae, we sequenced transcriptomes of tenebrionid genera from the Chilean Atacama Desert throughout. In addition to the transcriptomes of the Tenebrionidae from the Atacama Desert, the transcriptomes of a larger number of tenebrionid genera from other regions of the world were sequenced to improve taxon sampling for our transcriptome analyses. As a result, our dataset includes seven of the 11 described subfamilies and 47 tribes. We used these data to obtain the deduced amino acid sequences from 35 neuropeptide precursors per species. The suitability of neuropeptide precursor sequences for phylogenetic inferences was previously demonstrated in a proof-of-concept study (Bläser, Misof & Predel, 2020). This approach is relatively fast and simple as it is based on a limited set of easily identifiable and well conserved protein coding genes. In an alternate analysis using the same transcriptome dataset, the rather commonly used approach of compiling a large scale dataset of orthologous genes Figure 1 Overview of the study area in the Atacama Desert (shaded area). This region and the adjacent Andean Cordillera are home to about 34 genera of Tenebrionidae, whose phylogenetic relationships are analysed in this study. Also shown are selected representatives of individual genera. Number of Atacama species and total number of species within the genera are noted, respectively. The dotted blue line is the 4,000 m.a.s.l. contour line in the west and the dashed red line is the average annual rainfall isohyet of two mm. The lower panel shows an elevation profile within the study area, exemplified for a cross-section south of Antofagasta (green line) with tenebrionids typical of different elevation levels along this transect. Raster map made with Natural Earth.
Full-size DOI: 10.7717/peerj.14848/ fig-1 was performed. Both approaches, the concatenated dataset of neuropeptide precursors and the large scale dataset of orthologous genes were thus used in parallel to evaluate the relationships within the Atacama Tenebrionidae. These analyses resulted in maximum support for most, but not all branches, and enabled a first convincing assessment of the phylogenetic relationships of the Tenebrionidae of the Atacama Desert.

RNA extraction, cDNA library preparation and sequencing
Total RNA was extracted from samples stored in absolute ethanol or from individuals kept alive until tissue dissection. To avoid excessive RNA degradation in specimens stored in ethanol, head and pronotum of the beetles were separated from the rest of the body before transferring them into ethanol. In larger species, the body was additionally opened longitudinally with sterilized scissors. Without any treatment prior to storage in ethanol, the RNA was usually highly degraded, suggesting limited penetration of ethanol across the cuticle. Grinding of whole insects was avoided in order to enable the intestine to be removed later. Insects alive until tissue dissection were kept at 4 • C for 10 min before preparation. In most individuals (both ethanol and fresh material), after removal of the appendages (legs, elytra, antennae), the body was opened dorsally with sterilized scissors, the intestine was removed and the central nervous system (  Total RNA from each sample was quantified using Qubit RNA Assay Kit (Thermo Fisher Scientific) and subsequently subjected to quality control and RNA integrity number (RIN) as implemented in the Agilent 2100 Bioanalyzer system (Agilent Technologies, Waldbronn, Germany). Finally, RNA from CNS and remaining tissue from each sample were pooled together in equimolar concentrations for library preparations. This approach improved the detection of peptide precursor sequences, whose genes are mainly expressed in the CNS. Sequencing libraries (double-indexed) were prepared using 1 µg of total RNA with the Illumina R TruSeq R stranded RNA sample preparation kit (Cat.20020594; Illumina, San Diego, CA, U.S.A.). If the total RNA concentration was insufficient for standard library preparation, at least 2 ng of extract was pre-amplified using the Ovation RNA-Seq System V2 (NuGen, San Carlos, CA, USA). The library preparation of pre-amplified samples was performed according to the Nextera XT DNA sample preparation protocol (part no. 15031942 Rev. C). Subsequent sample preparation and sequencing was carried out at the Cologne Center for Genomics on an Illumina HiSeq 4000 and Illumina NovaSeq 6000 systems as described in Ragionieri & Predel (2020) with 75 bp or 100 bp paired end reads.

Transcriptome assembly, evaluation of cross-contaminations and statistics
Raw data (FASTQ files format) were filtered by removing adapter sequences and low quality using Trimmomatic 0.38 (Bolger, Lohse & Usadel, 2014 et al., 2011;Haas et al., 2013) with the read normalization option. All transcriptome assemblies were checked for potential cross-contaminations due to multiplex sequencing of several libraries using CroCo v1.1 (Simion et al., 2018) which removes potential sources of contamination using both transcriptome assemblies and the corresponding paired raw data (Table S1). This strategy uses sequence similarities and abundances to detect potential cross-contaminations. For closely related species that are analysed together, this can lead to an overestimation of cross-contamination (Simion et al., 2018). CroCo was run with the following settings: fold-threshold 2, minimum-coverage 0.1, overexp FLOAT 300, minimum percent identity between two transcripts to suspect across contamination 98%, minimum length of an alignment between two transcripts to suspect a cross contamination 180. Finally, we checked for and eliminated additional contamination of vector and linker/adapter using UniVec database (http://www.ncbi.nlm.nih.gov/tools/vecscreen/univec/). The transcriptomes assembled in this study lost on average about 2% of their sequence information due to the crosscontamination check. The quality and the completeness according to conserved singlecopy ortholog content of transcriptome assemblies were evaluated using the Perl script (TrinityStats.pl) included in Trinity and BUSCO v3 based on an Endopterygota obd9 dataset (Seppey, Manni & Zdobnov, 2019), respectively. The filtered transcriptome assemblies were submitted to NCBI Transcriptome Shotgun Assembly database (Table 1) and used for the large scale data set phylogenetic reconstruction.

Orthology assessment and alignment of neuropeptide precursors
Available amino acid sequences of neuropeptide precursors of Tr. castaneum and Te. molitor (Li et al., 2008;Veenstra, 2019;Marciniak, Pacholska-Bogalska & Ragionieri, 2022) were used as initial queries to search for orthologous sequences in the transcriptome assemblies. The assembled transcripts were analysed with the tblastn algorithms provided by NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi) or BioEdit version 7.0.5.3 (Hall, 1999). In case of missing data, precursor sequences of closely related taxa were used as alternative query sequences. Candidate nucleotide precursor gene sequences were translated into amino acid sequences using the ExPASy Translate tool (Artimo et al., 2012; http://web.expasy.org/translate/) with the standard genetic code. Orthologous neuropeptide precursor sequences were aligned using the MAFFT-L-INS-i algorithm (Katoh & Standley, 2013) (dvtditr (amino acid) Version 7.299b alg=A, model=BLOSUM62, 1.53, −0.00, −0.00, noshift, amax = 0.0); terminal sequences which were only found in few species were manually trimmed. The results were then manually checked for misaligned sequences using, e.g., N-termini of signal peptides and conserved amino-acid residues (cleavage signals, Cys as target for disulfide bridges) as anchor points. Individual amino acid alignments of each group of orthologous neuropeptide precursors were concatenated with catsequences 1.3 (https://zenodo.org/record/4409153#.YmJYT35Byot). The average evolutionary divergence for each neuropeptide precursor was calculated as in Bläser & Predel (2020). Briefly, overall mean distances (± standard error after 500 bootstrap generations) were computed with MEGA X (Kumar et al., 2018) implementing the Poisson correction model (Zuckerkandl & Pauling, 1965). Amino acid compositions and parsimony informative sites of the combined alignment were calculated using MEGA X.

Compilation of an orthologous gene dataset of Tenebrionidae
A Coleoptera orthologous reference gene set was compiled using OrthoDB v10. This approach provides reliable markers for phylogenomics (Misof et al., 2014;McKenna et al., 2019). Single copy genes shared across species of Coleoptera (Taxonomy ID: 7041) were selected for analysis. Orthograph (Petersen et al., 2017) was used to generate a profile hidden Markov model from the amino acid sequences of transcripts of each reference gene on the filtered transcriptome assemblies. Initially, we obtained 2,689 orthogroups (OGs) shared among Coleoptera, which were subsequently aligned using the MAFFT-L-INS-I algorithm (Katoh & Standley, 2013). Alignment ambiguities or spurious sequences in each OG were identified and removed using trimAL 1.2 (Capella-Gutierrez, Silla-Martinez & Gabaldon, 2009) with residue overlap threshold (-resoverlap 0.75) and sequence overlap threshold (-seqoverlap 90). With that approach, 947 out of 2,689 OGs were removed from the initial data set. Finally, all OGs were concatenated in a single partitioned super-alignment using catsequences.

Genome sequencing, assembly and identification of myosuppressin genes
Whole genome extraction was carried out using thoracic muscles of a single individual of Nycterinus abdominalis Eschscholtz, 1829 collected in Talcahuano

Phylogenetic analysis of neuropeptide precursors and a large scale orthologous gene dataset
FASTA files of aligned peptide precursor sequences were converted into PHYLIP and NEXUS formats using AliView 1.18-beta7 (Larsson, 2014). After defining the N-terminus of each neuropeptide precursor as starting partition, best-fit partitioning schemes and substitution models for subsequent phylogenetic analyses were predicted with ModelFinder (Chernomor, von Haeseler & Minh, 2016;Kalyaanamoorthy et al., 2017;Minh et al., 2021) implemented in IQ-TREE release 2.1.4b (Minh et al., 2020). Models and concatenated alignments for all analyses of both data sets are listed in Data S1 and S2. All phylogenetic analyses have been rooted using the Cleroidea Melyris sp. Bayesian inference (BI) analyses were run with MrBayes, with four runs, using eight chains and a sample frequency of 1,000 until convergence was achieved (PSFR value between 1.00-1.02) with a 10,000,000 generations (Ronquist et al., 2012). Maximum likelihood (ML) analyses were carried out using IQ-TREE 2.1.4b. ML analyses of both data sets were ran with the nearest-neighbour interchange search to consider all possible nearest-neighbour interchanges (-allnni) and branch support was evaluated with 1,000 ultra-fast bootstrap (UFBoot) (Hoang et al., 2017) and the Shimodaira-Hasegawa-like approximate likelihood ratio test (SH-aLRT) (Guindon et al., 2010). Trees were visualized using FigTree 1.4.2 (http://tree.bio.ed.ac.uk/) and designed in Inkscape 1.0 (https://inkscape.org/).
All analysed taxa of Alleculinae, Blaptinae, Diaperinae, and Stenochinae, as well as those taxa of Tenebrioninae that nest within the sister clade of Blaptinae (this clade is marked with an arrow in Fig. 2), have a distinct synapomorphy in common, namely an insertion of eight amino acids in the myosuppressin precursor ( Fig. 3A; see Data S3 for full sequences). This insertion does not result from differential transcription, but it is indeed manifested at the gene level. This could be verified by genome sequencing of a N. abdominalis specimen and a subsequent comparison of the myosuppressin gene structures (exons) of N. abdominalis and Tr. castaneum (Fig. 3B).

(B) Analysis of a large scale dataset of orthologous genes
The partitioned and concatenated alignment is composed of 1742 OGs with an overall length of 788,676 amino acid sites (Data S2). The best fitting models according to ModelFinder are listed for each partition in Data S2, which also contains the concatenated alignment. The topology of the resulting tree (Fig. 4) is largely congruent with that of the neuropeptide precursor data set. Differences are mainly observed for several of those branches with low support in the neuropeptide precursor tree (see Fig. S1): Salax as sister to a clade comprising Achanius, Arthroconus, Aspidolobus, Cordibates, Eremoecus, Hylithus, and Thinobatis; Praocis as sister to Falsopraocis + Physogasterini; Auladera + Nyctelia as sister to Callyntra + Psectrascelis; Sepidium as sister to Psammetichus + Eurychora; Alleculinae as sister to Scotobiini/Nycterinus/Zophobas + Scaurini; Nestorinus as sister to Heliofugus + Cuphotes Champion, 1887; and Isomira + Prionychus as sister to Omophlus + Heliotaurus. In addition, Discopleurus is sister to a main branch of Pimeliinae (Fig. 4), including, among other tribes, also the Stenosini; and Nycterinus changed its position and was recovered as sister to Diastoleus + Scotobius. In the large scale data set of orthologous genes, the branches with low support (Fig. S2) include that with Zophobas as sister to Nycterinus, Scotobius and Diastoleus (SH-aLRT = 8.2/UFBoot = 61). In both data sets, Corticeus has the same position, but the corresponding branches are very long.  Matthews et al. (2010), Bouchard et al. (2021 and Kamiński et al. (2020); Color coding for Tenebrionidae: Alleculinae, yellow; Blaptinae, grey; Diaperinae, light blue; Lagriinae, dark green; Pimeliinae, dark blue; Stenochinae, light green; Tenebrioninae, red. Atacama genera are marked with asterisks. The arrow marks the clade with a synapomorphy in the myosuppressin gene. Posterior probability (PP) and UFBoot (Bt) values are highlighted with circles on the nodes: black, above or equal to 0.95/95; grey, between 0.90−0.94/90-94; white, below 0.90/90. The detailed information on posterior probability/UFBoot values as well as the ML tree are provided in Fig. S1.

DISCUSSION
Transcriptomic information, mostly obtained from single individuals, was on the one hand used to obtain the amino acid sequences of 35 orthologous peptide precursors of genera of Tenebrionidae from the Atacama Desert and of selected taxa from other regions of the world. Due to their co-evolution with their corresponding receptors, neuropeptide sequences are particularly conserved and very well suited for a reconstruction of phylogenetic relationships at the intra-ordinal level (Bläser, Misof & Predel, 2020;Predel et al., 2012;Roth et al., 2009). Other advantages of using such datasets are the ease of ortholog assignment and the presence of unambiguous and highly conserved sequence motifs that facilitate a manual control of alignments generated by sequence alignment programs. The parallel analysis of the large scale dataset of orthologous genes revealed mostly the same topology as the neuropeptide precursor tree, with the exception of the few differences discussed below. The majority of Atacama genera cluster in three clades. Two of these clades belong to the subfamily Pimeliinae, which contains most of the desertadapted darkling beetles worldwide (Doyen, 1993;Kergoat et al., 2014b). In the Pimeliinae Pimelia/Akis were found to be the sister group to the remaining 17 analyzed tribes of Pimeliinae. The latter lineage consists of two clades with worldwide distribution, each containing a larger number of Atacama genera. One of these clades contains Nycteliini, Praociini and Physogasterini and forms a well-supported monophyletic group. This confirms previous morphological studies, which suggested Praociini, Physogasterini and Nycteliini as closely related taxa (Doyen, 1972;Doyen, 1993). These tribes are only known from arid regions of South America and are thought to be the sister group of North American Coniontini Waterhouse, 1858, Branchini LeConte, 1862and Asidini (Doyen, 1993. The Mediterranean Alphasida representing Asidini, was recovered in our analyses as sister to Praociini, Physogasterini and Nycteliini. Different from the most recent cladistic analysis of morphological characters in Nycteliini (Flores, 2000a), our analysis shows monophyletic Nycteliini as sister to Praociini + Phyogasterini. Within Nycteliini, which generally avoid the hyperarid core of the Atacama Desert, Pilobalia + Gyriosomus form the sister clade to the remaining Nycteliini. Physogasterini represent another very well supported monophyletic group in our analyses and include many species typical of the hyperarid core of the Atacama Desert. However, Praociini as defined in Flores (2000b) and Flores & Vidal (2009) appear polyphyletic with both data sets, this tribe requires a re-evaluation based on molecular data. From the four genera included here, Praocis and Falsopraocis are sister to Physogasterini, while Antofagapraocis and the central Chilean Gyrasida do not form a monophyletic group with Praocis and Falsopraocis. The genus Psammetichus, which is typical of hyperarid environments along the Coastal Cordillera of the Atacama Desert and the Pampa de Tamarugal, belongs to a sister clade of the above tribes. That clade also includes Sepidiini and Adelostomini, which do not occur in South America (Bouchard et al., 2021). Kamiński et al. (2022) suggested Sepidiini and Adelostomini as closely related tribes, considering the morphology of female terminalia and several genes. However, they did not place Elenophorini close to these tribes. Psammetichus was transferred to Elenophorini by Doyen & Lawrence (1979), a tribe that also includes Leptoderis (= Elenophorus Dejean 1821) of the western Mediterranean. Leptoderis was also included in our transcriptomic dataset, but the molecular data do not support an ancient link. In fact, Psammetichus was kept in Elenophorini in the past, although it was never found to be closely related to Leptoderis in Doyen's cladograms (Doyen, 1993). Also Ferrer (2015) doubts this relationship due to a number of morphological characters not shared between Leptoderis and the South American Elenophorini. In our tree, Leptoderis robustly nests within Stenosini, the latter represented by Chilean Discopleurus (within Stenosini only in the neuropeptide tree) and Hexagonochilus, and the palaearctic Dichillus as sister to Leptoderis.
The second major branch of Pimeliinae excl. Akis/Pimelia has Caenocrypticoides as sister to the rest. Caenocrypticoides is a well-established example of members of the same tribe (here Caenocrypticini;Endrödy-Younga, 1996) occurring in widely separated arid regions of Africa and South America, and thus probably representing a relict pattern that points to xerophilic ancestors before the break-up of Gondwana. The sister clade of Caenocrypticoides diverges into one lineage with diverse taxa having a wide distribution in the Palaeartic and Africa, but are not present in South America (Erodiini, Tentyriini, Zophosini, Adesmiini) and a second lineage with South American taxa belonging to Edrotini, Epitragini, Evaniosomini, Thinobatini, and Trilobocarini. The current placement of genera within these tribes is based on morphological characters (e.g., Doyen, 1993;Flores & Aballay, 2015). Although the exact position, particularly those of Arthroconus (Edrotini), Salax (Trilobocarini) and Achanius (Evaniosomini) could not be fully resolved with our data, it is obvious that none of the tribes is monophyletic. This South American clade was already mentioned by Doyen (1993) as a group ''not easy to fit with any classification'' using morphology and the classification at tribe level of the different genera have seen several changes over time (see e.g., Flores & Aballay, 2015). Doyen (1993) himself suggested transferring Achanius to the Edrotini ( =Eurymetopini Casey, 1907). The first split in this lineage separates Evaniosomus/Melaphorus/Aryenis (Evaniosomini) + Trilobocara (Trilobocarini) from the remaining taxa with maximum branch support. These remaining taxa include, among others, Achanius, Eremoecus, and Salax (Trilobocarini) and thus further genera of the aforementioned tribes and are separated in the neuropeptide tree into Geoborus/Nyctopetus (Epitragini) + Salax and a subclade which, in addition to Achanius, Arthroconus and Eremoecus, also includes Aspidolobus as another representative of the Epitragini. In the large scale dataset of orthologous genes, Salax is sister to all above mentioned taxa, including Geoborus + Nyctopetus. Finally, the well supported sister group relationship of Hylithus (Edrotini) and Thinobatis (Thinobatini) clearly argues against the supposed monophyly of Thinobatini which is only composed of the two genera included in our study (Doyen, 1993;Bouchard et al., 2021).
The sister group of Pimeliinae contains all other tenebrionid taxa analyzed in our study. The basal branching separates Lagriinae from the rest, which shows an early branching of Tenebrio + Bolitophagini and Tribolium + Melanimini. The remaining taxa split into the recently re-established Blaptinae (Kamiński et al., 2020) incl. Blapstinus from the Atacama Desert, and a diverse group of taxa including Stenochinae, Diaperinae, Alleculinae, and Tenebrioninae. Blapstinus appears to be the only tenebrionid genus from the Atacama Desert that has close relatives in North America. The corresponding subtribe Blapstinina Mulsant & Rey, 1853 is in fact restricted to Nearctic and Neotropical regions (Lumen et al., 2020;Kamiński et al., 2022). Monophyly of the analyzed taxa of Lagriinae, Blaptinae, Stenochinae, and Alleculinae was confirmed with maximum branch supports, respectively. On the other hand, polyphyly was evident for Diaperinae and Tenebrioninae (see also, e.g., Gunter et al., 2014;Kergoat et al., 2014b;Kamiński et al., 2020). Most taxa of the darkling beetles currently grouped in the subfamilies Alleculinae, Blaptinae, Diaperinae, Stenochinae, and Tenebrioninae have well-developed hindwings and do not show particular adaptations to hyperarid environments (Doyen, 1993). This does not apply to the Scotobiini, which represent the only endemic tribe of Tenebrioninae in arid South America (Matthews et al., 2010) and include the third cluster of tenebrionid genera in the Atacama Desert. In fact, three of the six genera of Scotobiini (Scotobius, Diastoleus, Ammophorus) inhabit the Atacama Desert and were included in our analysis. Within this clade Scotobius + Diastoleus is sister to Ammophorus in the neuropeptide tree, whereas in the large scale data set of orthologous genes Nycterinus replaces the position of Ammophorus. While the classification within Scotobiini of Diastoleus and the widespread Scotobius has been stable, the systematic position of the genus Ammophorus changed considerably over time. When Solier (1838) established the Scotobiini, he included Ammophorus in this tribe. Shortly afterwards Lacordaire (1859) transferred this genus to Nyctoporini Lacordaire, 1859 (Pimeliinae), where it remained for over 100 years (see, e.g., Kulzer, 1955;Peck, 2006;Peña, 1966b). Later, Vidal & Guerrero (2007) transferred Ammophorus to Elenophorini (Pimeliinae). Based on detailed analyses of morphological characters, Doyen (1993) and Silvestro, Giraldo-Mendoza & Flores (2015) proposed to return the genus to Scotobiini. The result of the neuropeptide tree fits the placement of Ammophorus within Scotobiini based on morphology (Silvestro, Giraldo-Mendoza & Flores, 2015). Also, they share a peculiar synapomorphy with the presence of dome-shaped placoid sensilla on the last segment of the antennae (Doyen, 1993). As sister of Scotobiini appears in the neuropeptide tree Zophobas Dejean, 1834 which is known only from Central and tropical South America (Ferrer, 2011). Nycterinus which is historically listed as the only South American genus within Amphidorini (see Doyen & Lawrence, 1979), belongs to the same monophyletic group in both data sets and was identified as sister to the above mentioned Scotobiini + Zophobas in the neuropeptide tree. Recent molecular phylogeny also showed Nycterinus as not belonging to the North American Amphidorini tribe, but rather to the South American Scotobiine clade which also includes Scotobiini and Zophobas (Johnston et al., 2022). The different results of the two data sets do not yet allow us to determine the specific position for Nycterinus.
The highly scattered appearance of the Tenebrioninae across the phylogenetic tree may question the reliability of our results. However, the topology does not show a mixture of taxa with poorly resolved sister group relationship, nor is it the result from particular poor taxon sampling. With the taxon-specific insertion of eight amino acids into the myosuppressin precursor (see Fig. 3) we have found a distinct synapomorphy at the molecular level clearly supporting Alleculinae, Blaptinae, Diaperinae, Stenochiinae, and a number of Tenebrioninae as a higher level monophyletic group. Based on morphological examinations, Doyen and Tschinkel speculated already in 1982 that Diaperinae, Stenochiinae, and Alleculinae could be derived offshoots of Tenebrioninae. In any case, it does not seem an easy task to redefine any clade as Tenebrioninae, except the one containing Tenebrio and Bolitophagini in our analyses.

CONCLUSIONS
Using newly generated transcriptome data, we were able to perform a comprehensive phylogenomic analysis of the tenebrionid fauna of the Atacama Desert and fill a gap in our knowledge of the highly diversified Tenebrionidae. The two datasets used for our analyses show only a few discrepancies that might be resolved by more extensive taxon sampling. The majority of Atacama genera are placed into three groups, two of which belong to typical South American lineages within the Pimeliinae. The suggested very close relationship of Psammetichus with the Mediterranean Leptoderis was not confirmed. Caenocrypticini including the Chilean Caenocrypticoides comprises a small group of genera present in southern Africa and (mostly) the Andean region of South America. These taxa display a combination of characters shared with various clades (Doyen, 1993). Our results provide the first evidence for a position of Caenocrypticoides as the sister of one of the main branches within Pimeliinae. While our data support the monophyly of the Nycteliini, Physogasterini and Scotobiini, this does not hold for the Atacama genera of Edrotini, Epitragini, Evaniosomini, Praociini, Thinobatini, Stenosini, and Trilobocarini. To clarify the relationships of these taxa, it is certainly useful to include more southern South American representatives in future analyses. In general, a detailed systematic revision of each of the latter groups appears necessary. As a side effect of our study, we have found a synapomorphy grouping Alleculinae, Blaptinae, Diaperinae, Stenochinae, and several taxa of Tenebrioninae, but not Tenebrio and Tribolium. This character, an insertion in the myosuppressin gene, defines a higher-level monophyletic group within the Tenebrionidae.